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ABSTRACT 

We discuss the properties of orbits within the influence sphere of a supermassive black hole (BH) , 
in the case that the surrounding star cluster is nonaxisymmetric. There are four major orbit families; 
one of these, the pyramid orbits, have the interesting property that they can approach arbitrarily 
closely to the BH. We derive the orbit-averaged equations of motion and show that in the limit of 
weak triaxiality, the pyramid orbits are integrable: the motion consists of a two-dimensional libration 
of the major axis of the orbit about the short axis of the triaxial figure, with eccentricity varying as 
a function of the two orientation angles, and reaching unity at the corners. Because pyramid orbits 
occupy the lowest angular momentum regions of phase space, they compete with coUisional loss cone 
repopulation and with resonant relaxation in supplying matter to BHs. General relativistic advance 
of the periapse dominates the precession for sufficiently eccentric orbits, and we show that relativity 
imposes an upper limit to the eccentricity: roughly the value at which the relativistic precession time is 
equal to the time for torques to change the angular momentum. We argue that this upper limit to the 
eccentricity should apply also to evolution driven by resonant relaxation, with potentially important 
consequences for the rate of extreme-mass-ratio inspirals in low-luminosity galaxies. In giant galaxies, 
we show that capture of stars on pyramid orbits can dominate the feeding of BHs, at least until such 
a time as the pyramid orbits are depleted; however this time can be of order a Hubble time. 
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1. INTRODUCTION 

Following the demonstration that self-consistent equi- 
libria could be constructe d for triaxial galaxy models 
()Schwarzschildl[l979l Il982| ) , observational evidence grad- 
ually accumulated for non-axisymme try on large (kilo- 
parsec) scales in ea,rly-type galaxies (iFranx et al.l 119911 : 
iStatler et all 120041: iCappellari et al.l [2007). On smaller 
scales, imaging of the centers of galaxies also revealed 
a wealth of features in the stellar distribution that 
are not consistent with axisymmetr y, including bars , 
bars-within-bars, and nuclear spirals phaw et al.l 119931 : 
lErwin fc Sparki[200l iSeth et al.l l2008'). In the^clei of 
low-luminosity galaxies, the non-axisymmetric features 
may be recent or recurring, associated with ongoing star 
formation; in luminous elliptical galaxies, central relax- 
ation times are so long that triaxiality, once present, 
could persist for the age of the universe. 

In a triaxial nucleus, torques from the stellar potential 
can induce gradual changes in the eccentricities of stel- 
lar orbits, allowing stars to find their way into the cen- 
tral BH. Gravitational two-body scattering also drives 
stars into the central BH, but only on a time scale of or- 
der the central relaxation time, which can be very long, 
particularly in the most luminous galaxies. Simple ar- 
guments suggest that the feeding of stars to the central 
BHs in many galaxies is likely to be dominated by large- 
scale torques rather than by two-body relaxation (e.g. 
iMerritt fc Poot][2001 . 

This paper discusses the character of orbits near a su- 
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permassive BH in a triaxial nucleus. The emphasis is 
on low-angular-momentum, or "centrophilic," orbits, the 
orbits that come closest to t he BH. Self-consistent mod- 
elling (iPoon fc Merrittll2004t ) reveals that a large fraction 
of the orbits in triaxial BH nuclei can be centrophilic. 

Within the BH influence sphere, orbits are nearly Ke- 
plerian, and the force from the distributed mass can 
be treated as a small perturbation which causes the or- 
bital elements (inclination, eccentricity) to change grad- 
ually with time. A standard way to deal with such 
motion is via orbit averaging (e.g. iSanders fc Verhuls^ 
119851 ). i.e., averaging the equations of motion over the 
short time scale associated with the unperturbed Ke- 
plerian motion. The result is a set of equations de- 
scribing the slow evolution of the remaining orbital el- 
ements due to the perturbing f o rces. This approach was 
followed by iSridhar fc Toumal ()1999D for motion in an 
ax ially-symmetric nucleus con taining a massive BH, and 
by iSambhus fc Sridharl ()2000D for motion in a constant- 
density triaxial nucleus. 

In their discussi o n of motion in triaxial nuclei, 
ISambhus fc Sridhail ()2000[ ) passed over one impor- 
tant class of orbit: the centrophilic orbits, i.e., or- 
bits that pass arbitrarily close to the BH. Exam- 
ples of centrop hiliic orbits include the t wo-di mensional 
"lens" orbits (jSridhar fc Toumal [19971 [l999h an d the 
three-dimensioiial "p yramids" (jMerritt fc Valluril Il999l: 
iPoon &: Merrittll2001|) . Centrophilic orbits are expected 
to dominate the supply of stars and stellar rem nants to 
a supermassive BH (e.g. IMerritt fc Poonll2004t ) and are 
the focus of the current paper. 

The paper is organized as follows. In ^we present a 
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model for the gravitational potential of a triaxial nuclear 
star cluster, which i s more general than that studied in 
iSambhus fc Sridhaii (|2000[ ). but which has many of the 
same dynamical features. Then in [J3]we write down the 
orbit-averaged equations of motion, and in Sj4l present a 
detailed analytical study of their solutions, with empha- 
sis on the case where the triaxiality is weak and the eccen- 
tricity large. In this limiting case, the averaged equations 
of motion turn out to be fully integrable. In ^ we derive 
the equations that describe the rate of capture of stars on 
pyramid orbits by the BH. Comparison of orbit-averaged 
treatment with real-space motion is made in ^ to test 
the applicability of the former. In [J7] we consider the 
effect of general relativity on the motion, which imposes 
an effective upper limit on the eccentricty. Sjg] discusses 
the connection with resonant relaxation: we argue that a 
similar upper limit to the eccentricity should character- 
ize orbital evolution in the case of resonant relaxation. 
Finally, in fJ5] we make some quantitative estimates of 
the importance of pyramid orbits for capture of stars in 
galactic nuclei. JfTOl sums up. 

2. MODEL FOR THE NUCLEAR STAR CLUSTER 

Consider a nucleus consisting of a BH, a spherical star 
cluster, and an additional triaxial component. An ex- 
pression for the gravitational potential that includes the 
three components is 



GM, ^ 



ro 



2-7 



(1) 



-2nGpt {T,x^+Tyy^ + T, 



The second term on the right hand side is the potential of 
a spherical star cluster with density p{r) = ps{r/ra) i"; 
the coefficient $s is given by 
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The scale radius tq may be chosen arbitrarily but it is 
convenient to set tq = rinfl, with rinfl the radius at which 
the enclosed stellar mass is twice M,: 
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The third term is the potential of a homogeneous triaxial 
ellipsoid of density pt ; this term can also be interpreted as 
a first approximation to the potential of a more general, 
inhomogeneous triaxial component. In the former case, 
the dimensionless coefficients (Tx,Ty,Tz) are expressible 
in terms of the axis ratios (jp, a) of the ellipsoid via el- 
liptic integrals ()Chandrasekhaii[l969() . The x{z) axes are 
assumed to be the long(short) axes of the triaxial figure; 
this implies T^ < Ty < T^. In what follows we will gener- 
ally assume pt <^ Psi'f'o), i-e. that the triaxial bulge has 
a low density compared with that of the spherical cusp 
at r = nnfl. 

3. ORBIT-AVERAGED EQUATIONS 

Within the BH influence sphere, orbits are nearly Ke- 
pleriarQ and the force from the distributed mass can be 

^ We consider general relativistic corrections in [[7] 



treated as a small perturbation which causes the elements 
of the orbit (inclination, eccentricity etc.) to change 
gradually wi th time. A stan dard way to deal with such 
motion (e.g. iSanders fc Verh ulst 1985) is to average the 
equations over the coordinate executing the most rapid 
variation, e.g., the radius. The result is a set of equations 
describing the slow evolution of the remaining variables 
due to the perturbing forces. 

We begin by transforming from Cartesian co- 
ordinates to action-angle variabl es in the K epler 
problem. Follow ing Sridhar & T oumal (|1999l ) and 
ISambhus fc Sridhaii (120001), we a dopt the Delaunay vari- 
ables (e.g. [Goldstein et alll2002[ ) to describe the unper- 
turbed motion. 

Let a be the semi-major axis of the Keplerian or- 
bit. The Delaunay action variables are the radial action 
/ — {GM,a)^''^, the angular momentum L, and the pro- 
jection of L onto the z axis Lz- The conjugate angle 
variables are the mean anomaly w, the argument of the 
periapse zu, and the longitude of the ascending node fl. 
In the Keplerian case, five of these are constants; the 
exception is w which increases linearly with time at a 
rate 



In terms of the new variables, the Hamiltonian is 



(3) 
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H = --f-^j +<^piI,L,Lz,w,uj,n); (4) 

the first term is the Keplerian contribution and $p, the 
"perturbing potential", contains the contributions from 
the spherical and triaxial components of the distributed 
mass. This transformation is completely general if we 
interpret the new variables as instantaneous (osculating) 
orbital elements. However if we assume that the per- 
turbing potential is small compared with the point-mass 
potential, the rates of change of these variables (again 
with the exception of w) will be small compared with 
the radial frequency i^^i and the new variables can be re- 
garded as approximate orbital elements that change little 
over a radial period P = 2-k /vj.. Accordingly, we average 
the Hamiltonian over the fast angle w: 
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The final term replaces the mean anomaly w by the ec- 
centric anomaly E, where r = a(l — ecosii^) and the 
eccentricity is e = -^/l — L? jP . After the averaging, 'K 
is independent of w, and / is conserved, as is the semi- 
major axis a. We are left with four variables and with 
$p as the effective Hamiltonian of the system. 
The spherically symmetric part of $p is 
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A good approximation to Fj{e) is 
-F-y(e) ~ 1 + ae 



2 a = ^i\^ ^: - 1 (7) 



¥r(4-7) 



which is exact for 7 = and 7 = 1; for < 7 < 2, 
< a < 3/2. When 7 > 1 and e is close to 1, a better 
approximation is 

VTT r(3-7) 

(8) 
We adopt the latter expression in what follows. Good 
approximations are 
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(9b) 



Similar expressions can be found in 

Ivanov. Polnarev fc Sahal ([2005) and 



Polyachenko. Polvachenko fc Shukhman (2007) . 

Expressions for the orbit-averaged triaxial harmonic 
potenti al (excluding the spherical component) are de- 
rived in iSambhus fc Sridhail (|2000f ) . Adopting their no- 
tation, the orbit-averaged potential in our case becomes 
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$p = $s — (1 + a - a'r) + 2TTGptT^ a^ x (10a) 
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(lOd) 



The shorthand Sx,Cx has been used for sin a;, cos a;. We 
have defined the dimensionless variables £ = L/I and 
£z = Lz/I , both of which vary from to 1; the orbital 
inclination i is given by cosz = £z/£ and the eccentricity 
by e2 ^ l-£2. 

The first term in equation (jlObp , which arises from the 
spherically-symmetric cusp, does not depend on the an- 
gular variables, and has the same dependence on £'^ as the 
corresponding term in the harmonic triaxial potential. 
So we can sum up the coefficients at £'^ and renormalize 
the triaxial coefficients €b.c to obtain the same functional 
form of the Hamiltonian as in the purely harmonic case. 
Dropping an unnecessary constant term (depending only 
on a) and defining a dimensionless time r = Vpt, where 
Vp is characteristic rate of precession. 
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(12) 



we obtain the dimensionless Hamiltonian and the equa- 



tions of motion describing the perturbed motion: 
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dr dw ' dr d£ ^ dr dil ' dr 9£ 

The renormalized triaxiality coefficients are e^ 
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If there were no spherical component (^4 = 0), this 
would reduce to the purely harmonic triaxial case studied 
by Sambhus & Sridhar (2000)0 Adding the spherically 
symmetric cusp increases the rate of periapse precession, 
while at the same time reducing the relative amplitude of 
the triaxial terms; otherwise the form of the Hamiltonian 
is essentially unchanged. 

A more transparent expression for the precession fre- 
quency i^p is 
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where M{a) denotes the mass enclosed within radius 
r = a. From equation (I13bp . the precession rate of an 
orbit in the spherical cluster is v^ = \dw/dt\ = Mvp = 
3\/l — e2 i/p. Near the BH influence radius, it is clear 
that i/p sa j/^; hence the orbit-averaged treatment, which 
assumes only one "fast" variable, is likely to break down 
at this radius. 

We note that i/^ — > as e — >■ 1. For the very eccentric 
orbits that are the focus of this paper, the rate of preces- 
sion is much lower than for a typical, non-eccentric orbit 
of the same energy. This will turn out to be important, 
since the slow precession allows torques from the triaxial 
part of the potential to build up. 

4. ORBITAL STRUCTURE OF THE MODEL POTENTIAL 
4.1. General remarks 



The orbit-averaged Hamiltonian (|13a|) describes a dy- 
namical system of two degrees of freedom. The trajec- 
tories must be obtain ed by numerical integration of the 
equations of motion (|13bp . We begin by making some 
qualitative points about the nature of the solutions. 

In the absence of the triaxial terms in equation (|13ap , 
the effect of the distributed mass is to rotate the pe- 
riapse angle -07 in a fixed plane; this steadily rotating 
elliptic orbit fills an annulus. The addition of a weak tri- 
axial perturbation changes the rate of in-plane precession 
slightly, and also causes the orbital plane itself t o cha nge, 
as described by the last two terms in equation (J13b|) . 

In general, two angular variables zu and fl can either 
librate around fixed points or circulate, giving rise to 
four basic families of orbits (Figure [T|). Solutions to the 
equations of motion that are characterized by circula- 
tion in both zu and fl correspond to tube orbits about 
the short axis (SAT). Motion that circulates in vj but 
librates in il corresponds to tube orbits about the long 



^ Equation (12) of lSambhus fc Sridharl lj200Ci 1 for £ lacks a minus 
sign in front of the first term. 
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Fig. 1. — Four classes of orbits around a BH in a triaxial nucleus. Left column: time dependence of the dimensionless angular momentum 
I (top/blue) and its component iz along the short axis of the figure (bottom/red). Middle column: argument of the periapse zu. Right 
column: angle of nodes Q. 



- V„,Vn 
1= 



^max— 'y 3 ^ 










to th e tube orbits that a re generic to the triaxial geom- 
etry ()Schwarzschildlll979[ ) . A subclass of the SAT orbits 
corres ponds to motion that circulates in fl and librate s 
in w (jSambhus fc Sridhari [2?M IPoon fc MerrittJl2nOH ). 
These orbits resemble cones, or saucers; similar orbits ex- 
ist also at r ^ r j nfi in oblate or nearly oblat e potentials 



0.01 0.1 1 

Fig. 2. — Dependence of the characteristic frequencies Ut^ (in- 
plane precession) and i/q (nodal precession) on the value of the 
dimensionless angular momentum £. Top (open) symbols: v^; bot- 
tom, (filled) sym,bols: uq. Magenta boxes: LATs; red triangles: 
SATs; yellow diamonds: saucers; green circles: pyramids. For 
large £ u^ cc (. and vq^ <^ Ut^, but for sufficiently low £ these two 
are comparable, which gives birth to the pyramid a nd sa ucer orbits. 
Vertical line denotes the threshold in £ (equa tion I36I I. horizontal 
line the characteristic frequency u^o (equation |22} . Triaxiality co- 
efficients were set to ec = 10"'^, ej, = 0.4ec- 

axis (LAT). Both types of orbit are qualitatively similar 



(|Richstonelll982tlLees fc Schwarzschildlll992D . 

If the degree of triaxiality is small {eb,c "C 1), then 
as noted above, the dominant effect of the distributed 
mass is simply to induce a periapse shift, at a rate 
i/ro = dw/dr = —2>l. If the additional mass is much 
less than the mass of the BH, then on short time scales 
(comparable to the radial period v~^) the orbit resem- 
bles a nearly closed ellipse. On intermediate time scales 
(of order the precession time y^^)^ a steadily-rotating 
elliptic orbit fills an annulus in a fixed plane. On still 
longer time scales Vq^ , the orbital plane itself changes 
due to the torques from the triaxial potential. Similar 
considerati ons give rise to the concep t of vector resonant 
relaxation (jRauch fc Tremainill996l) . 

The foregoing description is valid as long as the angu- 
lar momentum is not too low. Since the precession rate 
is proportional to €, for sufficiently low £ the intermedi- 



Orbital structure of triaxial nuclei 



ate and long time scales become comparable (Figure [5]). 
As a result, the triaxial torques can produce substantial 
changes in i (i.e. the ecc entric ity) on a precession time 
scale via the first term in (jf 3b|) , and the circulation in w 
can change to libration. This is the origin of the pyra- 
mid orbits, which are u nique to the triaxial geometry 
(|Merritt fcValluril [19991) . 

4.2. Pyramid orbits 

Of the four orbit families discussed above, the first 
three were treated, in the orbit-averaged approximation, 
by Sambhus & Sridhar (2000). The fourth class of orbits, 
the pyramids, are three-dimensional a nalogs of the two- 
dimen sional "lens" orbits discussed by lSridhar &: Toumal 
(|1997f l. also in the context of the orbit-averaged equa- 
tions. An important property of the p yramid orbits is 
that i can come arbitrarily close t o zero ()Poon fc Merritti 
I200H) and lMerritt fc PoonI (J2004). This makes the pyra- 
mids natural candidates for providing matter to BHs at 
the centers of galaxies. 

Pyramid orbits can be treated analytically if the fol- 
lowing two additional aproximations are made: (1) the 
angular momentum is assumed to be small, £^ ^ 1; (2) 
the triaxial component of the potential is assumed to 
be small compared with the spherical component, i.e. 
CbiCc <C 1. As shown below, these two conditions are 
consistent, in the sense that ^^ax ~ ^b,c for pyramid or- 
bits. 

Removing the second-order terms in ej,, ec and in i'^ 
from the orbit-averaged Hamiltonian ()13a[) . we find 

3 5 

H = --f + - [ec(l - cf)s% + eb{c^sn + CiSr^cn)'^] 

(15) 
where again Ci = cosi = ^z/^- 

Because an orbit described by ([T5|) is essentially a pre- 
cessing rod, one expects the important variables to be 
the two that describe the orientation of the rod, and 
its eccentricity. This argument led us to search for ex- 
act solutions to the equations of motion in terms of the 
Laplace-Runge-Lenz vector, or its dimensionless counter- 
part, the eccentricity vector, which point in the direction 
of orbital periapse. 

We therefore introduced new variables Cx, Sy and Cz'- 



of equations (J16p and using equations (J13b[) , we find 

Bx = 3^(sinn7 cos ft + cosro sin 51 cosi), (18a) 
By — 3^(sinn7 sin 51 — costu cos 51 cosi) (18b) 

where Cx = dex/dr etc. Taking second time derivatives, 
the variables describing the orientation and eccentricty 
of the orbit drop out, as desired, and the equations of 
motion for Cx and By can be expressed purely in terms of 
Cx and Cy: 

ex^~ex6iH + 3f) (19a) 

= -e, [30ec -QH~ aOe^el ~ 30(ec - eb)el], 

ey = -ey6{H + 3e^~^eb) (19b) 

^-ey [30ec ~6H^15eb~ 30e,el - 3G(e, - eb)el] . 

From equations ()18|) . {ex,ey) = implies I — Q, i.e. the 
eccentricity reaches one at the "corners" of the orbit. 
These define the base of the pyramid. Defining (ea:Oieyo) 
to be the values of (e^:, Cj,) when this occurs, vthe Hamil- 
tonian has numerical value 
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(20) 



Equations ([T9| have the form of coupled, nonlinear os- 
cillators. Given solutions to these equations, the time 
dependence of the additional variables (^,^2,tx7, 51) fol- 
lows immediately from equations (1161) and (jl8|) : 



^x ' ^y y^x^y ^x^y ) 
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9(1 -el -el) 9ie. + ej, + ej, 

tz y^x^y ^x^y } / '^^ 
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y^x^x ~r ^y^y ) 
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and a quite lengthy expression for 51 which we choose not 
to reproduce here. 

In the limit of small amplitudes, {ex,ey) ^ 1 (which 
corresponds to 7J ~ f^c), the oscillations are harmonic 
and uncoupled, with dimensionless frequencies 



Ct ^cos-cu cos 51 



By = sin w cos i cos 51 
ez = sintA7 sini 



sin nj cos i sin 51 (16a) 

costu sin 51 (16b) 

(16c) 



which correspond to components of a unit vector in the 
direction of the eccentricity vector. Of these, only two 
are independent, since el + e^ + el = 1. In terms of these 
variables, the Hamiltonian (|15|) takes on a particularly 
simple form: 



H 



€c^^ 



(cc - eb)el 



(17) 



and Cy, which describe the orientation of the 



As expected, the Hamiltonian depends on only three vari- 
ables: 6; 
orbit's major axis, and the eccentricity i 

To find the equations of motion, we must switch to a 
Lagrangian formalism. Taking the first time derivatives 



'^xo = VlScc , i^yQ = \/l5(ec - ef,). 



(22) 



The apoapse traces out a 2d Lissajous figure in the plane 
perpendicular to the short axis of t he triaxial figure. Thi s 
is the base of the pyramid (e.g. iMerritt fc Valluril ()1999() , 
Figure 11). The solutions in this limiting case are 

ex{T)^exaCOs{iyxQT + (l)x), (23a) 

ey{T) = eyQCOs{i^yoT + (j)y), (23b) 

^^iT)=^lo si^^i^xoT + 0x) + ^lo sin^ii'yor + (t>y) 

where (f'xji'y ^-re arbitraty constants and 

^xO — VxQexol'i, f-ya — VyVlS-yo/'i- (24) 

Figure [3^ plots an example. 

Equations (1231) describe integrable motion. Remark- 
ably, it turns out that the more general (anharmonic, 
coupled) equations of motion ([T9)) are integrable as well. 
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Fig. 3. — A pyramid orbit, in three approximations. Each orbit has the same {exO,£yo) = (0.5,0.35). (a) The simple harmonic oscillator 
(SHO) approximation, equations I I23I I. valid for small l'^, {ei,,ec) and (e2;o,eyo)- (b) From equation II19I I. which does not assume small 
(63:0,6^0). (c) From the full orbit- averaged equations H13| l. which does not assume small £^, e or {ej:0,eyo)- Aside from the fact that 
the latter orbit is fairly close to a 5 : 2 resonance, the correspondence between the physically important properties of the approximate 
orbits is good. The triaxiality parameters are (e;,, Ec) = (0.0578, 0.168), corresponding to a pyramid orbit with a = O.lro in a nucleus with 
triaxial a xis r atios (0.5,0.75), density ratio ptiro)/ps{ro) = 0.1, and 7 = 1. The frequencies for the SHO case are UxO = 1.5 9, Uyo = 1.28 
(equation [22J ; frequencies for planar orbits with the same e^ and Sy amplitudes are 1.48 and 1.24 respectively (equation I27I I. 



The first integral is H; an equivalent, but nonnegative, 
integral is U where 



U = 15ea-6H 



^oe^ + (e'+e^ + e')- (25) 



The second integral is obtained after multiplying the first 
of equations (fT9|) by 15ecex, the second by 15(ec — ti,)ey, 
and adding them to obtain a complete differential. The 
integral W is then 



W^vlM + 



2 2 



vine 



xo^i) + i^loiel + vie. 



V V 



2 4 

K,r,e. 



n Z 2 

lU + vlr 



2 2 



u 



•'yo- 



yO'^y) 

(26a) 
(26b) 



The existence of two integrals {U, W), for a system with 
two degrees of freedom, demonstrates regularity of the 
motion. 

Regular motion can always be expressed in terms of 
action-angle variables. The period of the motion, in each 
degree of freedom, is then given simply by the time for 
the corresponding angle variable to increase by 27r. We 
were unable to derive analytic expressions for the action- 
angle variables corresponding to the two-dimensional mo- 
tion described by equations ([19)). However the periods 



of oscillation of the planar orbits (e^ = or ej, = 0) 
described by these equations are easily shown to be 



VxoP{exi^)^-^K 



^xOl 



(ej 



VyQP{eyo)=AK{el^i) (e 



0), 
0) 



(27) 



where K{a) is the complete elliptic integral: 



K{a) 



7r/2 



1 



>-l/2 



dx. 



For small a, K ^ 7r/2 and P ~ 2ti/vq. As a — )■ 1, 
K -^ oo] this corresponds to a pyramid that precesses 
from the z axis all the way to the (x, y) plane. The 
oscillator is "soft": increasing the amplitude increases 
also the period. Figure [3] shows comparison of orbits with 
the same initial conditions, calculated in three different 
approximations . 

Pyramid orbits can be see n as analogs of regul ar box 
orbits in triaxial potentials (jSchwarzschildl [l979[ ). with 
three independent oscillations in each coordinate. Like 
box orbits, they do not conserve the magnitude of the 
sign of the angular momentum about any axis. The dif- 
ference is that a BH in the center serves as a kind of "re- 
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Fig. 4. — Regions in the U — W plane occupied by t he d ifferent 
orbit famihes. Lower (da rk b lue) line (a), equation II29II; upper 
(red) curve (d), equati on l|34|l : green line (c), eq uati on l|30|l : light 
blue line (b), equation II31I I: red points, equation I I35I I: blue points, 
equation g2]l. Plotted for e^ = 0.002, ec = 0.008. 

fleeting boundary", so that a pyramid orbit is reflected 
by 180° near periapsis, instead of continuing its way to 
the other side oi x — y plane as a box orbit would do. 

4.3. The complete phase space of eccentric orbits 

While our focus is on the pyramid orbits, the low- 
angular-momentum Hamiltonian ()f 5|) also supports or- 
bits from other families. In this section, we complete the 
discussion of the phase space described by equation ([T5|) . 
by delineating the regions in the U — W plane that are 
occupied by each of the four orbit families (Figure [J]) . 

Pyramids and LATs both resemble distorted rectan- 
gles in the ex,ey plane. The corner points of this region 



correspond to e^ 



0. Evaluating the two integrals 



at a corner (denoted by the subscript 0) gives 

U = via elo + vIq el„ + e^ g , (28a) 

{U-el,)elo- (28b) 



W = vto e^o + ^to 4o 



The difference between pyramids and LATs arises from 
the last term: corner points of pyramid orbits correspond 
to £2 = and hence (from (|21a|) ) to e^,o = 0. For LATs 
the condition is, conversely, e^ g > 0, and Cz.o — (hence 
CyQ = 1 — e^p). Analyzing these expressions, we find that 
for pyramids and LATs the lower and upper boundaries 
for W given U are 



W^iy'yoU, 



(29) 
(30) 



and the boundary between pyramids and LATs is given 
by 

W^{vl, + vl,)U-i^l^vl^. (31) 

Pyramids lie above and to the left of this line in the 
U — W plane, while LATs are below and to the right. 
The intersection of this line with (p9)) and (|30)) occurs at 
the points 



U = vIq, W = Vyo, 



(32a) 

(32b) 

These points constitute the leftmost bound for LATs and 
the rightmost bound for pyramids respectively. 



Short-axis tubes and saucers resemble distorted rect- 
angular regions in the ex^&z plane. Again, the corner 
points (with subscript 0) are defined to have e^ = 62 = 



and Cy = 0, with e^ > 0, and therefore 



V- 



2 1 



-yQ ■ 



W^vloU+ {i^loelo + vIq - ul^) e. 



-yo- 



(33a) 
(33b) 



Both these families have W > eJJ , i.e. lie above the line 
((30|) . SAT orbits intersect the plane Cz = 0, so we can 
set CxQ = 1 in ([33|l . (Alternatively, for SATs, both angles 
circulate, so we can set w — U, = Q, which again gives 
ex = l)- We then find that SATs lie below the line ([5T]) . 
On the other hand, saucers never reach e^ = (since 
for them sin^ w > 0), so that they lie above the line (|3T|) . 
To obtain the upper limit for W at fixed U , we substitute 
e^g from the first equation in ([55)1 in the second, and 
then seek a maximum of W with respect to Cxo at fixed 
U . This gives 

W^yl^U+{U + i^l^-vl^f/A. (34) 



(35a) 
(35b) 



This curve intersects (|30l) and (|3T|) in the points 
U^i^lo + ^vo. W^iy^oU + i^to, 



which define the left- and rightmost bounds for the saucer 
region. 

All these criteria are summarized in Figure |4l In par- 
ticular, pyramid orbits exist in the following cases: 

• for < if < |e6 they coexist with LATs; 



• for < if < |(ec 



Cb) they coexist with SAT 



above these values they are the only population for 
H < ^Cr, which is the maximum allowed value of 



H. 



2'=cj 



• below H < pyramids do not exist (this is easily 
seen from equation (|15p : since the term in square 
brackets is always non-negative, it is impossible to 
have 1^ = when ii < 0) . 

Figure \5\ shows Poincare surfaces of section for fl = 
n/2 and < ii < |ef,. The three families of orbits are 
delineated. 

Since H is an integral of the motion, the maximum 
allowed value of i'^ can not exceed 



AH) = 



-^J-^^«i(5e.-2ii). (36) 
S + 4ec - £6 3 



The latter approximate expression is immediately seen 
from the simplified Hamiltonian (jlSp . while the former 
comes from the exact Hamiltonian (J13p . However, it does 
not follow that an orbit with a sufficiently low instanta- 
neous value of the angular momentum is necessarily a 
pyramid: both tube families can also have arbitrarily 
low £. The principal distinction is that any pyramid or- 
bit can achieve arbitrarily low £ (that is, the lower bound 
is £ = 0), while tube orbits always have < £^i„ < £^ 
(however small fmin may be, it is strictly positive). 
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Fig. 5. — Poincare section for i^vo plotted at SI = 7r/2 for en- 
ergy H = 0.02 (efc = 0.99~2 - l.^c = 0.96"^ - 1) showing the 
three possible types of orbit: LATs, pyramids and SATs/saucers. 
This figure adopts the s ame potential parameters as Figure 5b in 
ISambhus fc Sridh"arl 120001 ). but those authors chose H = —0.02 
which precludes pyramid orbits. Boundaries are marked by the 
same letters as in Figure HI 



We now return from the simplified Hamiltonian (jlSp to 
the fuh Hamiltonian (J13p . i.e. we no longer require e to 
be small. The full Hamiltonian retains all the qualitative 
properties of the simplified system but requires numerical 
integration of the equations of motion ()13|) to determine 
orbit classes. 

To quantify the overall fraction of pyramid orbits in a 
given potential, one should uniformly sample the phase 
space for all four variables and determine the orbit class 
for each initial condition. From equation p6|) . we can 
restrict ourselves to values of £'^ < ^niax(O) = 3+4^"-^ 
(but we must take care not to filter out initial conditions 
corresponding to iJ < 0). 

We calculated the proportions of the i^^g^^-restricted 
fraction of phase space occupied by each family of orbits. 
Initial conditions were drawn randomly for 10'* points 
(with uniform distribution in i^ G [^--^max]: '^^ ^z G [0--^]j 
and in -07, ri G [0, §])• The proportions were found to 
depend very weakly on Ec if fc ^ 1. To elucidate the 
dependence on eb/ec (the degree of triaxiality) we took 
15 values in the range (0.001 - 0.999). 

We found that the relative fraction rj of pyramids 
among low-£ orbits is almost independent of ec (Figure|6]): 



0.28 




(37) 



The fraction of pyramids among all orbits is fj = 
^^max(O) = f^c??- For comparison, the left panel of 
Figure [6] shows the results obtained using the simpli- 
fied Hamiltonian ([15]) and the analytical classification 
scheme described above, while the middle panel, made 
for Ec = 0.1, shows almost the same behavior, with the 
addition of a small number of chaotic orbits. We note 
that for Ec ~ 1 the phase space becomes largely chaotic. 

These estimates of the relative fraction of pyramid or- 
bits are directly applicable to a galaxy with an isotropic 
distribution of stars at any energy. This assumption may 
not be valid, for example, in the case of induced tangen- 
tial anisotropy following the m erger of supermassive BHs 
(jMerritt fc Milosavlievi3l2005D . 

One can ask a different question: if we know the in- 
stantaneous value of an orbit's eccentricity and orienta- 



tion, what can we conclude about the orbit class? It is 
clear that without knowledge of the derivatives of ex,y 
the answer will only be probabilistic. It turns out that 
the probability p for an orbit with "sufficiently high" ec- 
centricity (i.e. with £'^ < ^,^,iax) to be a pyramid depends 
mostly on the z component of the eccentricity vector: 

0.7w44^(l — ^) e\-^ (here the normalization comes 



p«U.^^4-(l ^ 

from the total number of pyramids among low-i orbits). 
That is, an orbit lying in the plane defined by the long 
and intermediate axes of the potential is certainly not a 
pyramid, and the highest probabilitiy occurs for orbits 
directed toward the short axis. 

4.4. Large i limit 

In the previous sections we considered the case et^c ^ 1 
and i'^ ~ Ec, which allowed a simplification leading to 
integrable equations. 

In the opposite case, when et^c ^ ^^ ^ 1, the frequency 
of in-plane precession, v^, is much greater than the rates 
change of 17 and i (Figure [21) . In this limit we can carry 
out a second averaging of the Hamiltonian (|13ap . this 
time over vj. Thus 



1 r 3 



(38) 



+ ^^ [e, (4 + cfcf^ + E, (5 - 3f) (1 - c^)] . 

On timescales T > v^^ the orbit resembles an annulus 
that lies in the plane defined by the angles i and ft. The 
only remaining equations of motion are those that de- 
scribe the change in orientation of the orbital plane: 



§— |(5-3£^)(l"C,^)^2n, 
^§ = |(5-3£^)4c,-|(5-3£^ 



(39) 



One expects the natural variables in this case to be the 
components of the angular momentum: 



£x =^sin j sinil, 
£y =£sin JcosJ7, 
£z—£ cos i 



(40) 



and £l + £y + £l — £^ — constant. In terms of these 



^ ' 2 4^2 



variables, the Hamiltonian is 

[e,{£^-£l)+e4£'-£l)]. 

(41) 
After some algebra, one finds the equations of motion: 

(42) 



4 = -^(6c-E.)(5-3^2)^ 



g2 ' 



4 = 1 (5 -3£^)^, 
4 = -|(5-3f2^44 



£2 



(only two of which are independent). These can be writ- 
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Fig. 6. — Proportions of the restricted part of phase space (defined by fi < ^maxC)) that are occupied by the major orbit families: LATs, 
pyra mids, SATs, and SAT/saucers, as a function of t he ra tio €i,/ec- Left: analytic estimates from the simplified orbit-averaged Hamiltonian 
Us} for ec — > 0; middle: orbit-averaged Hamiltonian ||13|I for ec = 0.1; right: real-space integration for orbits with semimajor axis a = ri^fj 
(equal to the BH influence radius) and ec = 0.1. 



ten 



d7 



:TX^, 



5-3^2 



(43) 



2P 




Conservation of the Hamiltonian (PT|) implies 
f-bl^l + fc^z = constant = C. 

This is an elliptic cylinder; the axis is parallel to the ix- 
axis, and the ellipse is elongated in the direction of the 
£j,-axis. In addition, we know that 

(l+(-l+il^ constant = f 

which is a sphere. So, the motion lies on the intersec- 
tion of a sphere with an elliptic cylinder. There are two 
possibilities. 

1. P > Cjth- In this case, the cylinder intersects the 
sphere in a deformed ring that circles the £a;-axis. This 
corresponds to a LAT orbit. 

2. (? < C/eb- In this case, the locus of intersection is a 
deformed ring about the f^-axis. This orbit is a SAT. 
In other words, precession of the angular momentum vec- 
tor can be either about the short or long (not intermedi- 
ate) axes of the triaxial ellipsoid. 

5. CAPTURE OF PYRAMID ORBITS BY THE BH 

As we have seen, pyramid orbits can attain arbitrar- 
ily low values of the dimensionless angular momentum 
i. The BH tidally disrupts or captures stars with angu- 
lar momentum less than a certain critical value L,, or - 
in dimensionless variables ~ i, = L,/I{a). We can ex- 
press (.t in terms of the capture radius rt , the radius at 
which a star is either tidally disrupted or swallowed. For 
BH masses greater than ^ 10^ Mq, main sequence stars 
avoid disruption and rt ~ rschw = IGM^ji?; for smaller 
M,, tidal disruption occurs outside the Schwarzschild ra- 
dius; e.g. at the center of the Milky Way, rt ~ lOrschw 
for solar-type stars. Defining rt = Orschw and writing 
L\ K, GM,rt, then gives 



e 



y'Schw 

a 



iQ-^e 



/ M. 



V lO^M© / V 1 Pc 



(44) 



We note the following property of the pyramid orbits: 
as long as the frequencies of e^ and Cy oscillation are in- 
commensurate, the vector {ex, Sy) fills densely the whole 
available area, which has the form of distorted rectangle. 
The corner points correspond to zero angular momen- 
tum, and the "drainage area" is similar to four holes in 
the corners of a billiard table. 

Unless otherwise noted, in this section we adopt the 
simple harmonic oscillator (SHO) approximation to the 
{^XT^y) motion, that is, we use the simplified Hamil- 
tonian (fT5l) and its solutions (|23|) : these orbits have 



e^ -I- e^ <C 1 and they form a rectangle in the Cx — Cy 
plane, with sides 2ex, 2ej,. As long as the motion is inte- 
grable, the results for arbitrary pyramids with ex,ey ^ 1 
will be qualitatively similar. Quantitative results may be 
obtained by numerical analysis and are presented near 
the end of this section. 

Figure [7] shows a two-torus describing oscillations in 
{cx, Cy) for a pyramid orbit. In the SHO approximation, 
solutions are given by (|23|) . If the two frequencies Uxo, Vyo 
are incommensurate, the motion will fill the torus. In 
this case, we are free to shift the time coordinate so as 
to make both phase angles {4>i,4'2) zero, yielding 

^^ {t) = ^lo sin^ [vxot) + ilo sin^ (lyyor) (45a) 



= 4oSm^^i+^^ 



yO^ 



(45b) 



where 9i = Vxqt, O2 — i^yoT- (In the case of exact com- 
mensurability, i.e. miVxo + rn2Vyo = with (7711,7712) 
integers, the trajectory will avoid certain regions of the 
torus and such a shift may not be possible.) In the 
SHO approximation, Vxo — vT5e^, VyQ = •\/l5(ec — £&) 
((22|) . More generally, integrable motion will still be rep- 
resentable as uniform motion on the torus but the fre- 
quencies and the relations between £ and the angles will 
be different. 

Stars are lost when (-{61,02) < t,- Consider the loss 
region centered at {01,62) = (0,0). This is one of four 
such regions, of equal size and shape, that correspond to 
the four corners of the base of the pyramid. For small 
it, the loss region is approximately an ellipse. 






^y0n2 < 1 
^^2 < 1. 



(46) 
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Fig. 7. — Two-torus describing oscillations of [e^^e-y) for a pyra- 
mid orbit. The ellipses correspond to regions near the four corners 
of the pyramid's base where (.<(.,. In the orbit-averaged approx- 
imation, trajectories proceed smoothly along lines parallel to the 
solid lines, with slope tana = Vy/u^. In reality, successive periapse 
passages occur at discrete intervals, once per radial period. 



The area enclosed by this "loss ellipse" is 

D2 



l^^t 



yO 



(47) 



There are four such regions on the torus; together, they 
constitute a fraction 



1 
/i= - 



TT ^lO^yO 



(48) 



of the torus. 

Stars move in the (6*1,6*2) plane along lines with slope 



tana = Vyo/vxo, at an angular rate of \ v'^q + i^yQ- Since 

periapse passages occur only once per radial period, a 
star will move a finite step in the phase plane between 
encounters with the BH. The dimensionless time between 
successive periapse passages is At — 2'KVp/i'r. The angle 
traversed during this time is 



A6' ^2-K(vj,lvr)^vlQ 



'yQ- 



(49) 



The rate at which stars move into one the four loss el- 
lipses is given roughly by the number of stars that lie 
an angular distance A0 from one side of a loss ellipse, 
divided by At. 

This is not quite correct however, since a star may 
precess past the loss ellipse before it has had time to 
reach periapse. We carry out a more exact calculation 
by assuming that the torus is uniformly populated at 
some initial time, with unit total number of stars. To 
simplify the calculation, we transform to a new phase 



plane defined by 



v.: 



/, ,2 f2 I , .2 tf2 
Y '^xO^xO ^ yO yO 

— VynC-xof^ydf^l + t^xoixoiyQ(^2 



(50) 

(51) 



With this transformation, the phase velocity becomes 



^ = (^'o4o + ^yo4o) 



1/2 



?9 = 



and the loss regions become circles of radius £,. 
angular displacement in one radial period is 



At/^ = 2^(z.p/z.,) {i^l^ilo 






(52) 
The 

(53) 



The density of stars is {4:Tr'^ixoiyo)'~^- 

At any point in the (V', "&) plane, stars have a range 
of radial phases. Assuming that the initial distribution 
satisfies Jeans' theorem, stars far from the loss regions 
are uniformly distributed in x where 



Vr 



(54) 



here P = lix jv^ is the radial period, r^ is the periapse 
distance and Vr is the radial velocity. The integral is 
performed along the orbit, hence x ranges between and 
1 as r varies from Vp to apoapse and back to r^. [x — 
w mod 27r, where w is mean anomaly). 

Figure [8] shows how stars move in the (x, ip) plane 
at fixed -d. The loss region extends in tfj a, distance 
2y/£l — z?2, from -0in to i/'out- Stars are lost to the BH if 
they reach periapse while in this region. 

Two regimes must be considered, depending on 
whether Aip is less than or greater than ^out — ''/'in- 

1. Ai/) < ■(/'out — V'in (Figure [8^). In one radial period, 
stars in the orange region are lost. One-half of this region 
lies within the loss ellipse; these are stars with i < io but 
which have not yet attained periapse. The persistence of 
stars inside the "loss cone" is similar to what occurs in 
the case of diffusional loss cone repopulation, where there 
is also a "boundary layer" , the width of which depends 
on t he ratio of the relaxatio n time to the radial period 
(e.g. ICohn fc KulsrudI (|1978f )). The other one-half con- 
sists of stars that have not yet entered the loss region. 
The area of the orange region is equal to the area of a 
rectangle of unit height and width Ai/; since stars are 
distributed uniformly on the (x, ip) plane, the number of 
stars lost per radial period is equal to the total number 
of stars, of any radial phase, contained within Aip. 

2. Alp > V'out — V'in (Figure [SId). In this case, some 
stars manage to cross the loss region without being cap- 
tured. The area of the orange region is equal to that of a 
rectangle of unit height and width V'out — "i/in • The num- 
ber of stars lost per radial period is therefore equal to 
the number of stars, of arbitrary radial phase, contained 
within f/out - "0111 = 2a/^2 _ ^2_ 

To compute the total loss rate, we integrate the loss 
per radial period over ■&. It is convenient to express the 
results in terms of q where 



q^ 



A^P 
24 



xO^xO 



4- ;y2 / 



yO^yO- 



(55) 
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Fig. 8. — Trajectories of stars in the (V'lX) plane as they en- 
counter a loss region from left to right, defined as ipi^ < V" 5: '/'out- 
X increases from at periapse, to 1/2 at apoapse, to 1 at subsequent 
periapse. Trajectories are indicated by dashed lines. Stars are lost 
if they reach periapse while inside the loss region. Stars within the 
orange region are lost in one radial period, (a) A^ < ^jj, — i/^out; 

(b) Alp > i>in - V'out. 



q -^ 1 corresponds to an "empty loss cone" and q ^ 1 
to a "full loss cone". However we note that - for any 
q < 1 ~ there are values of 1} such that the width of the 
loss region, V'out — "0111, is less than Aip. In terms of the 
integral W defined above ((28)) . q becomes simply 



q 



PVp 



W. 



(56) 



Unlike the case of coUisional loss cone refilling, where 
q = q{E) is only a function of energy, here q is also a 
function of a second integral W . Pyramid orbits with 
small opening angles will have small W and small q. 

The area on the (f/', "&) plane that is lost, in one radial 
period, into one of the four loss regions is 



AV-di?- 



2 / (V'out 



Vin)dl9 



(57) 



where 



is the value of •& where A?/; = V'out 



(58) 
V'in; for Q > 1, 




Fig. 9. — Illustrating the area of the (1/), "d) phase plane that is 
lost into the BH each radial period. The circle centered at (0, 0) is 
the loss region corresponding to one corner of the pyramid orbit; its 
radius is I, . Regions marked in bold denote the area of the torus 
that is lost in one radial period, for 9 < 1 (Ai/)i) and g > 1 (A^2). 
While the number of stars lost per radial period is proportional to 
the marked areas, the region on the torus from which those stars 
come is more complicated since it depends also on an orbit's radial 
phase (Figure (8)1. 

i?c = 0. For g < 1, the area integral becomes 



Aql, / dd + A ^/ei - d^dd 

r-l 



-Aqlly/l-q^ + Ail / dx^Jl - x^ 

-^(lU + 2qyJ\-q^ - 2 arcsin yJl-qA 
■Milfiq), 



/(<?) = 2^1 -g^ + Y arcsin((7) 



(59) 



and for q > 1 it is tt^^. The function f(q) varies from 
/(0) = 1 to /(I) = 7r/4 « 0.785. 

The area on the phase plane that is lost each radial 
period can be interpreted in a very simple way geomet- 
rically, as shown in Figure [SI 

Considering that there are four loss regions, the in- 
stantaneous total loss rate J-, in dimensionless units, is 



•^ - /(g) ^2 ^ 'f , V'^'o ^lo + ^lo ^lo 



M 4q/(q) 



for < g < 1, 



Pl/p TT 



(60a) 



T = q- 



,,2 />2 , 2 f2 _ 



27r4o4oV -0 -°" y°^y° 2^2 4g£ 



PVr, 



for q > 1. 



yQ Vp 

(60b) 



The second expression for the loss rate, equation (|60b| . 
can be called the "full-loss-cone" loss rate, since it corre- 
sponds to completely filling and empyting the loss regions 
in each radial step (Figure [5]). Note that the loss rate for 
q < 1 \s ^ q times the full-loss-cone loss rate. A similar 
relation holds in the case of coUisionally repopulated loss 
cones (Cohn fc Ku lsrud 1978) . 

The inverse of the loss rate J-" gives an estimate of 
the time idrain required to drain an orbit, or equivalently 
the time for a single star, of unknown initial phase, to 
go into the BH. In this approximation, the loss rate re- 
mains constant until t = idrain at which time the torus 
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Fig. 10. — Proportion of phase space (£^ < ^maxC^)) that is occu- 
pied by the major orbit families: LATs, pyramids, SATs, saucers, 
and chaotic orbits, as a function of semimajor axis a based on 
real-space integrations. Triaxiality p aram eters are ej, = 0.5ec and 
Cc = 0.12/[0.2-|-(a/rinfl)~^] (equation [TTJ, corresponding to a den- 
sity cusp with 7 = 1 an d £c = 0-1 at a = rjng . For a > rinfj most 
\ow-l orbits are chaotic IjPoon fc Merritt|[200!l 'l . 

is completely empty. In reality the draining time will 
always be longer than this, since after ~ 1 precessional 
periods, some parts of the torus that are entering the loss 
regions will be empty and the loss rate will drop below 
equation (|60|) . For A?/' > V'ln — V'out, the downstream 
density in Figure [H integrated over radial phase, is eas- 
ily shown to be 1 — q'^^^Jl — -ff^ j i\ times the upstream 
density while for A?/; < V'ln ~ V'out the downstream den- 
sity is zero. Integrated over i?, the downstream depletion 
factor becomes 

1_ ^_ 713^(1 + ^) + _Lsin-iyT^^ (61) 

for g < 1 and 1 — 7r/4g for q > 1; it is for g = 0, 
^ 0.215 for g = 1 and 1 for g — > oo. For small g, the 
torus will become striated, containing strips of nearly- 
zero density interlaced with undepleted regions; the loss 
rate will exhibit discontinous jumps whenever a depleted 
region encounters a new loss ellipse and the time to to- 
tally empty the torus will depend in a complicated way 
on the frequency ratio v^^jvy and on ^,. For large g, the 
loss rate will drop more smoothly with time, roughly as 
an exponential law with time constant ~ tdrainO 

We postpone a more complete discussion of loss cones 
in the triaxial geometry to a future paper. Here we make 
a few remarks about pyramids with arbitrary opening 
angles, i.e. for which e^^Oj ej/O are not required to be small. 

For each orbit one can compute ^, the fraction of the 
torus occupied by the loss cone (equation H5|) . by nu- 
merically integrating the equations of motion (1131) and 
analyzing the probability distribution for instantaneous 
values oi P-: V(P- < X) <y^ X - €^j„, where £^j„ allows 
for a nonzero lower bound on £^. Almost all pyramids 
have imin = 0, but some of them happen to be resonances 
(commensurable v^, and Vy) and hence avoid approach- 
ing ^ = 0. This linear character of the distribution of 
€^ near its minimum corresponds to a linear probability 
distribution of periapse radii (7^(?'peri < r) ex r), which 
is natural to expect if we combine a quadratic distribu- 
tion of impact parameters at infinity with gravitational 

^ This was the approximation adopted by Merritt & Poon (2004). 



focusing (see equation 7 of iMerritt fc PoonI ()2004[ )V 

The coefficient [i for each orbit is calculated as 'P{f' < 
£l). As seen from equation (|48)) . the smaller the extent 
of a pyramid in any direction, the greater fi - this is 
true even for orbits with large Cxo or Cyo- While fi varies 
greatly from orbit to orbit, its overall distribution over 
the entire ensemble of pyramid orbits follows a power 
law: 



Vf^ifi > r) « {Y/^i„^^ny 



2?7' 



(62) 



Vfi is the probability of having /i greater than a certain 
value and rj is the fraction of pyramids among all orbits 
(|37p . The average /z for all pyramid orbits is therefore 
JL = 2fj,min, and the average fraction of time that a ran- 
dom orbit of any type and any £ spends inside the loss 
cone is Jifi ~ l^ (almost independent of the potential pa- 
rameters efc and Ec) ~ the same number that would result 
from an isotropic distribution of orbits in a spherically- 
symmetric potential. 

6. COMPARISON WITH REAL-SPACE INTEGRATIONS 

We tested the applicability of the orbit-averaged ap- 
proach by comp aring the orbit-averaged equations of mo- 
tion, equation (J13b|) , with real-space integrations of or- 
bits having the same initial conditions (and arbitrary ra- 
dial phases). The agreement was found to be fairly good 
for orbits with semimajor axes a < O.lrjnfl: about 90% 
of the orbits were found to belong to the same orbital 
class, and the correspondence between values of fj, and 
^min ^^^ ^^^'^ quite good for individual orbits. Averaged 
over the ensemble, the proportion of phase space occu- 
pied by the different orbital families, as well as the net 
flux of pyramids into the BH, is almost the same for the 
two methods. However, at larger radii, the relative frac- 
tion of pyramids and saucers decreases (Figure |6] (right), 
nop . Since the maximum possible angular momentum for 
orbits with a given semimajor axis a grows faster than 
^/GM^a, this means that the fraction of pyramid orbits 
among all (not just low-^) orbits is even smaller. For 
orbits with semimajor axis a > rinfl the frequency of ra- 
dial oscillation becomes comparable to the frequencies of 
precession, and when these overlap, orbits tend to be- 
come chaotic. (Weakly chaotic behavior starts earlier). 
So low-^ orbits with a > l.Srinfi are mostly chaotic, as 
seen from Figure [TUl confirming that regular pyramid or- 
bits (along with saucers) exist only within BH sphere of 
influence (|Poon fc Merrittl[200l . 

7. EFFECTS OF GENERAL RELATIVITY 

In the previous sections we considered the BH as a 
Newtonian point mass. In general relativity (GR), the 
gravitational field of the BH is more complicated, and 
this will affect the behavior of orbits with distances of 
closest approach that are comparable to Vg = GM,/c^ . 

For a non-spinning BH, the lowest order post- 
Newtonian effect is advance of the periapse, which acts in 
the opposite sense to the precession due to an extended 

* We note that saucer orbits also exist in potentials with high 
central concentration of mass, such as logarithmic potential studied 
in e.g. Lees fc Schwarzschild! II1992I ). 
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mass distribution. The GR periapse advance is 

67r GM, 



An7 = 



c? (1 — e-^)a 



(63) 



per r adial period, with c the speed of light (jWeinbergj 
|1972|) , making the orbit-averaged precession frequency 



3GAf. 



(64) 



We can approximate the effects of this precession by 
adding an extra term to the orbit-averaged Hamiltonian 



_ VGRf."^ 3GM, Vr rschw M, 



c^a Vp a M{a) 



(65a) 
(65b) 



This is equivalent to adding the term >c/P to the equa- 
tion of motion for tn, i.e. to the right hand side of 
dvo/dr = dH/dL When (. = i'crit, where 



(I) 



1/3 



(66) 



the precession due to GR exactly cancels the preces- 
sion due to the spherical component of the distributed 
mass. Since the angular momentum of a pyramid or- 
bit approaches arbitrarily close to zero in the absence of 
GR, there will always come a time when its precession 
is dominated by the effects of GR, no matter how small 
the value of the dimensionless coefficient x. 

We again restrict consideration to the simplified Hamil- 
tonian (|T5|) . vahd for eb,c <C 1, ^^ < Cc, now with the 
added term due to GR. This Hamiltonian may be rewrit- 
ten as 



H- 



2 

P{i) 



T^ce, 



+ ni^c- e&)e. 



Vl^xj Gj/j : 



(67) 



where P and Q denote the expressions in the first and 
second sets of square brackets. The minimum of P{£) 

occurs at € = ^crit^ 



Prr 



^81>f2/8. 



(68) 



The function Q can vary from to some maximum value 
Qmax due to the limitation that e^ -I- e^ < 1 . 
Two differences from the Newtonian case are apparent. 

1) For each value of {cx, ey) (and therefore Q), there are 
now two allowed values of i. One of these is smaller than 
^crit while the other is greater f Figure [TT]) . 

2) Both the minimum and maximum values of £ - both of 
which correspond to the maximum value of P (Figure[TT|) 



are attained when Q = 0, i.e. when e^ 



0. 



The maximum of Q corresponds to i = IcHt- In the 
Newtonian case, the minimum of £ corresponds to the 
maximum of Q. 

7.1. Planar orbits 

We first consider orbits confined to the y — z plane 
{Cx = 0, hence $7 = TT/2^iz = throughout the evolu- 
tion). Namely, we start an orbit from w = it/2 {ey = 0) 



0.03 



0.02 



0.01 




0.15 



Fig. 11. — Illustrating the allowed variations in angular momen- 
tum (. for orbits in the presence of general relativistic preces sion . 
The solid (red) curve represents the function P{t) (equation I67II : 
the dashed (blue) parabola is the same function in the Newtonian 
case {k = 0). If >r 7^ 0, P(l) has a minimum at ^^rit (equation I66II . 
Orbits make excursions along the curve P{t.) in the range from a 
certain value Pmax to Pmax — Qmax (here Qmax is given for the case 
of planar orbits of ^ 17. II and equals i{^c — ^b) = ^yo/*^)' ^^ during 
such an excursion £ does not cross ^criti then the orbit resides on 
one branch of P{i), typically remaining regular. Otherwise it flip s 
to the other branch, reaching lower values of ^niin (equation I79|l . 
becoming a (typically chaotic) LAT or pyramid orbit. 

and £ = Iq. In the absence of GR, such an orbit would 
be a LAT for £0 > |eti and a pyramid otherwise. 
The Hamiltonian and the equations of motion are 



r- 



'^-|^»+|-5'^+?+^»'-<'=») 



6 



sin2ti7 , 



~3£- 



(70) 



The orbit in the course of its evolution may or may not 
attain zu — (modTr). If it does, then the angle zu circu- 
lates monotonically, with w ^ 0. In Figure [11] the con- 
dition tjj = corresponds to reaching the lowest point in 
the P{£) curve, £ = 4rit- Whether this happens depends 
on the value of £q : since the orbit starts from Q = and 
P ~ P{£o), it can "descend" the P{£) curve at most by 
Qmax — Vy[)l^- If this condition is consistent with reach- 
ing P(^crit), the orbit will flip to the other branch of the 
P{£) curve. The condition for this to happen is 



t-o± 



i;o± 



r)^crit 



"vo. 



«crit 



(71) 



£0+ and £q- are the upper and lower positive roots of 
this equation. 

If £q > £q+, the orbit behaves like a Newtonian LAT 
(Figure [121 case a): it has w < and £ > ^crit- If 
£0 < £0-, the orbit is again a LAT, but now it precesses 
in the opposite direction (tij > 0) due to the dominance 
of GR, and £ never climbs above ^crit (Figure fT2l case e). 
In these cases the condition w = gives the extremum 
of £, which is found from equation (j69p : 



3^2 1 ^ 


-2^° + 4 


r,^oxtr,L 1 /, 

^ tlcxtr.L 



6 



(72) 



This extremum appears to be a minimum (^cxtr.L < ^0) 
if £0 > £0^ and a maximum (£extr,L > -^o) if ^0 < ^o-- 
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Pyramid orbits are those that reach ^ = fcrit- '!^ 
changes sign exactly at €ciit , but the angular momentum 
continues to decrease beyond the point of turnaround, 
reaching its minimum value only when w returns again 
to 7r/2, i.e the z-axis. The two semiperiods of oscillation 
are not equal: the first {£ > ^crit and rij < 0) is slower, 
the other is more abrupt (Figure [121 cases &, d). In 
effect, the orbit is "reflected" by striking the GR angular 
momentum barrier. After the orbit precesses past the 
z axis in the oposite sense, the angular momentum be- 
gins to increase again, reaching its original value after 
the precession in w has gone a full cycle and the orbit 
has returned to the z axis from the other side. 

If £o — ^crit, there is no oscillation at all - the GR 
and extended mass precession balance each other exactly 
(Figure [HI case c) . For £o £ ^crit the orbit precesses in 
the opposite sense to the Newtonian precession. 

We can find the extreme values of £ by setting ^ = in 
equation ([70]). This occurs for -co — tt/2, i.e. for Q = 
or P{£) = P{£o). This gives 



t-oxtr,P 



(^^l + 8£i,J£l-iy (73) 



If £o > £cTit, this root corresponds to the minimum £, 
with £q the maximum value; in the opposite case they 
exchange places. For >c ^ 3£q this additional root is 



2/^ 



2 K 

3B 



(74) 



Thus the minimum angular momentum attained by a 
pyramid orbit in the presence of GR is approximately 
proportional to k. Note the counter-intuitive result that 
the pyramid orbit with the widest base (largest £o) comes 
closest to the BH. 

Figure [T3| shows the dependence of the maximum and 
minimum values of £ on £q for the various orbit families. 

7.2. Three-dimensional pyramids 

In the case of pyramid orbits that are not restricted to 
a principal plane, numerical solution of the equations of 
motion derived from the Hamiltonian (|65ap are observed 
to be generally chaotic, increasingly so as >c is increased 
(Figure [H]) . This may be attributed to the "scattering" 
effect of the GR term k/1 in the Hamiltonian, which 
causes the vector {ex,ey) to be deflected by an almost 
random angle whenever £ approaches zero. In the limit 
that the motion is fully chaotic, H remains the only inte- 
gral of the motion. The following argument suggests that 
the minimum value of the angular momentum attained 
in this case should be the same as in equation ()73p . 

Suppose that the Hamiltonian (|571) is the only integral 
that remains. Then the vector {ex,ey) can lie anywhere 
inside an ellipse 

5 

Q{ea:,ey) = - [ecel + (cc - eb)ey] < Qmax , (75) 

whose boundary is given by 
5 

Qmax = XCc - ff - -Pmin- (76) 

^ T. Alexander has suggested that these be called "windshield- 
wiper orbits." 
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Fig. 12. — Planar y — z orbits, solutions of equation lj69|l in 
a potential with ec = 10~^,ei, = ec/2,x = 10^*, started with 
ro = tt/2 and different ^q: (a) 0.104, (b) 0.103, (c) 0.0322, (d) 
0.006, (e) 0.0058. The first two orbits lie close to the separatrix 
between LATs and pyramids, £o+ = 0.10392 I I71II : the third is 
the stationary orbit with Iq = ^crit; and the last two lie near the 
separatrix between pyramids and GR-precession-dominated LATs, 
^0— = 0.005845. Top panel shows the evolution of ^(r), bottom 
panel shows tjj{t). For pyramid orbits (b-d), the angle zu librates 
around it/2, and £ crosses the critical value £crit; tube orbits (a,e) 
have 07 monotonically circulating, and £ is always above or below 



This ellipse defines the base of the "pyramid" (which 
now rather resembles a cone). As in the planar case, 
the maximum and minimum values of £ are attained not 
on the boundary of this ellipse (i.e. the corners in the 
Newtonian case), but at Cx = Cy = 0, where Q = and 
P attains its maximum. These values are given by the 



roots of the equation P{£) 

M^ - (5ec - 2H)£ 



H, or 

6£l^0. 



(77) 



The two positive roots of this cubic equation are given 
by 



L-inin.inax 



= 3^^ 



/TT 

2H sin - 
V6 

9:x 



±1 



(78) 



r''^^°\(5..-2i7)3/2 



The plus sign in the argument of the sine function gives 
^max while the minus sign gives ^min- These two values 
are linked by a simple relation: 



(v/l + 8 



■'^^-max 



/3-1 



O -C™. 



(79) 
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Fig. 13. — Minimum and maximum values of l for a series of 
orbits with initial conditions f = £oi f^ = f/S, ^z = 0, n = 7r/2. 
Potential parameters are tc = 10~^,ei, = ec/2,x = 10~*. The 
straight li ne i s £ = £o; dashed line is the extremum for pyra mids , 
equation II73I I . These two curves intersect at £(.rit (equation I66II . 
where they exchange roles. For i > £o+ a^nd i < io— feauation l71l l 
the orbit is a tube, and the minimum (or maximum) is given by 
equation l|72|l . Dotted grey line shows the leading frequency of 
zu oscillations, u^ X 10^^; for high-^ orbits i^^ ~ 3i, for orbits 
dominated by GR precession t^ro sa 2k /£^, Letters denote the 
position of orbits shown in Figure [121 

where the latter approximate equality holds for x <^ 
^max- I'^ tti^ same approximation 



2m: 



5ec - 2H 



max 



5€r - 2H 



(80) 



Equation (|73p for planar pyramids is a special case of 
this relation where (,q — ^max- 

The ellipse (|75l) serves as a "reflection boundary" for 
trajectories that come below i « fcrit- If this happens, 
the vector (e^cgy) is observed to be quickly "scattered" 
by an almost random angle (Figure [HI right, denoted by 
the red segments) , similar to the rapid change in w that 
occurs in the planar case (Figure [T2|). Roughly speaking, 
all pyramid orbits and some tube orbits (those that may 
attain (. < i?crit) will be chaotic. 

The distinction between pyramids and chaotic tubes is 
in the radius of this ellipse: pyramids by definition have 
a fixed sign of e^, or e^ + e^ < 1, which means that the 
ellipse ([75]) should not touch the circle e^ + e^ = 1. Hence 
pyramids have Q < |(ec — Cfc), and 



5 5 



Pn- 



(81) 



This condition is different from the one described in Sj4l 
even in the case >c = = Pmin, since now pyramids do 
not coexist with LAT orbits. 

The condition for LATs to be chaotic (i.e. to pass 
through 4rit) is Pmin + Qmax < fcc - i?. For LATs the 
ellipse ([75)1 always intersects the unit circle, so this con- 



dition can be satisfied for —P„ 



< H < ^€b- Pn 



However, this is a necessary but not sufficient condition 

® A small fraction of the "flipping" orbits, especially those that 
oscilla te n ear f^rit (close to the lowest point on the P{i) curve of 
Figure [TTJ, may retain regularity by virtue of being resonant. 



for a chaotic LAT: some orbits from this range do not at- 
tain £ < f crit because of the existence of another integral 
of motion besides H (that is, they are regular). 

Finally, we consider the character of the motion when 
the precession is dominated by GR, as would be the case 
very near the BH. This is equivalent to staying on the 
left branch of P(£), with dn.max < ^o- < 4rit dZD)- In 
this limit there is a second short time scale in addition 
to the radial period, the time for GR precession. This 
situation is similar to the high-i' case (S I4.4p . in the sense 
that we can carry out a second averaging over va and 
obtain the equations that describe the precession of an 
annulus due to the triaxial torques. The orbits in this 
case are again short- or long-axis tubes. The only differ- 
ence from ^ 14.41 is that we have to add the term —x/i to 
the averaged Hamiltonian (|4T|) , but since £ is constant in 
this approximation, the equations of motion for £z,^ do 
not change. These very-low-£ regular tube orbits can be 
easily captured by the BH, however their number is very 
small and we do not consider them when computing the 
total capture rate. 

We argue in §8 that the conservation of £ for orbits in 
this limit can have important consequences for resonant 
relaxation. 

7.3. Capture of orbits by the BH in the case of GR 

The inclusion of general relativistic precession has the 
effect of limiting the maximum eccentricity achievable by 
a pyramid orbit. However, if fmin < £,, an orbit can still 
come close enough to the BH to be disrupted or captured. 



Introducing the quantity w 



1 



2k 



w ■ 



t» €, 06c 



> 



2H 



, we can write 

3 Vr £» 



(82) 



where we used (|64[ I65bp and set i/ = as a lower limit 
for pyramid orbits (orbits with the smallest H have the 
largest ^max and the smallest ^min)- Comparison with 
equation ([55|. with W < (15ec)^, shows that 



w ' 



Stt 



9 



(83) 



Roughly speaking, the condition that stars be captured 
(w < 1) is equivalent to the statement that the loss cone 
is full {q > 1). This is not a simple coincidence: a full 
loss cone implies that for low-£ orbits the mean change 
of £ during one radial period (~ I'o^o t^p/i^r) is of order 
£,, while the condition w = 1 requires that for the lowest 
allowable £, the GR precession rate (l64|) is comparable to 
the radial frequency. These two conditions are roughly 
equivalent. 

We can express this necessary condition for capture 
in terms of more physically relevant quantities. Writing 
equation p4ap as 



i^p ^ I M{a) 
Vr ^ 2 M, ' 

and approximating rjufl of equation ^ as 

GM, _ 



(84) 



(85) 



with a the one-dimensional stellar velocity dispersion at 
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Fig. 14. — Three pyramid orbits with the same initial conditions (£ = 0.05,^2 = 0.02, sj = Q = n/2;ec = 0.01, ej, = 0.005) and three 
values of the GR coefficient k (equation [6501 . Left: x = (regular); middle: x = 10 ~^ (weakly chaotic); right: x = 10~^ (strongly 
chaotic). The green ellipse marks the maximal extent of the {ex,ey) vector, equation Il75|l . i.e. £ = ^crit, equation l|66| l: red segments 
correspond to i < <?criti blue to £ > ^^rit a^nd to the nonrelativistic case. 



Tinfl, the condition w < 1 becomes 



1 



/Qc 



7-7/2 



<1. 



(86) 



The Milky Way BH constitutes one extreme of the BH 
mass distribution. Writing 9 « 10 (solar-mass main 
sequence stars), a ~ 10^ km s^^, and 7 = 3/2 gives 






(87) 



e.g. for a = O.lro « 0.3 pc, Cc ^ 10"'^ is required for stars 
to be captured. This is a reasonable degree of triaxiality 
for the Galactic center. 

At the other mass extreme, we consider the galaxy 
M87, for which a « 350 km s'^ and 9^3. Setting 
7 = 0.5, corresponding to a low-density core, we find 



.1/3. 



> 10- 



A dimensionless triaxiality of order unity is reasonable 
for a giant elliptical galaxy. 

In § [3] we estimate the loss rate using the expression 
(j60bp for the full loss cone rate, with the modification 
that the fraction of time ^ an orbit spends inside the 
loss cone is now given not by (liSl) . but by 



e 



fi'. 



(•2 _ f2 ■ 
•^0 <-min 



(89) 



This relation implies that the instantaneous value of i'^ 
is distributed uniformly in the range [^min--^max]- This 
is a good approximation for chaotic pyramid orbits and 
a reasonable (within a factor of few) approximation for 
chaotic LATs (and also for regular orbits). 

In the next section we point out the importance of 
the angular momentum limit for the rate of gravitational 
wave events due to inspiral of compact stellar remnants. 

8. CONNECTION WITH "RESONANT RELAXATION" 

Resonant relaxation (RR) is a phenomenon that arises 
in stella r systems exhibitiiig certain regularities in the 
motion ([Ranch fc Tremainill996HHopman fc Alexandeii 
l2006al) . Due to the discreteness of the stellar distribu- 
tion, torques acting on a test star from all other stars 



do not cancel exactly, and there is a residual torque that 
produces a change in the angular momentum: 



dL 
'dt 



N = L, 

a 



Nm, 



M, 



-2ttP- 



(90) 



(here m is the stellar mass, P 
period, Lc 



2t:/ fir is the radial 



\/GM,a is the angular momentum of a cir- 
cular orbit with radius a, and N is roughly the number 
of stars within a sphere of radius a). In a non- resonant 
system this net torque changes the direction randomly af- 
ter each radial period, but in the case of near-Keplerian 
motion, for example, orbits remain almost the same for 
many radial periods, so the change of angular momentum 
produced by this torque continues in the same direction 
for a much longer time, the so-called coherence time tcoh, 
until the orientation of either the test star's orbit or the 
other stars' orbits change significantly. If this decoher- 
ence is due to precession of stars in their mean field, then 



ic 



^M 



M, P 

m N 



(91) 



where the relevant precession time is that for an orbit of 
average eccentricity. 
The total change of L during tcoii is 



(AL) 



coh 



Gm 



■tr 



N 



(92) 



On timescales longer than icoh the angular momentum 
experiences a random walk with step size (AL)coh and 
time step icoh- The relaxation time is defined as the time 
required for an orbit to change its angular momentum by 
Lc, and hence it is given by 



tRB.,s « ( ^ 



t. 



coh 



P^. (93) 



The above argument describes "scalar" resonant relax- 
ation, in which both the magnitude and direction of L 
can change. On longer timescales, precessing orbits fill 
annuli, which also exert mutual torques; however, since 
these torques are perpendicular to L, they may change 
only the direction, not the magnitude of L. This effect is 
dubbed vector resonant relaxation (VRR) , and its coher- 
ence time is given by the time required for orbital planes 
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to change. In a spherically symmetric system the only 
mechanism that changes orbital planes is the relaxation 
itself. Q Hence for VRR the coherence time is given by 
setting |dL/di| — Lc/tco\i in equation (|90|) : 



icoh = io,VRR 



M, P 



m 



M 



and the relaxation time, equation ()93p . becomes 

Af. 

iflfl,u ~ ^n,VRR ~ P ■ 



m\/N 



(94) 



(95) 



which is ^ V iV times shorter than the scalar relaxation 
time. 

We begin by comparing RR timescales with timescales 
for orbital change due to a triaxial background potential. 
Consider a star on a (regular) pyramid orbit confined to 
the x—z plane. It experiences periodic changes of angular 
momentum i = L/Lc with frequency < Vxoi^p (|22p and 
amplitude £xo S ^xo/'i ([24l) . Hence, the typical rate of 
change of angular momentum is 



'^-^ ~ r 2 /q 



Lr 5er. 



Nm 

m7 



27rP- 



(96) 



Comparison with ([M|) shows that the rate of change of 
angular momentum due to unbalanced torques from the 
other stars (RR) is greater than the rate of regular pre- 
cession if evW < 1. However, the coherence time for 
RR is a typical precession time of stars in the cluster, 
Vp^, whereas pyramids change angular momentum on a 

longer timescale (i/pVTSe)"^. On the other hand, in the 
case of RR the angular momentum continues to change 
in a random- walk manner on timescale longer than icoh, 
while in the case of precession in triaxial potential its 
variation is bounded. 

Next we consider VRR, which corresponds to changes 
in orbital planes defined by the angles fl and i = 
aiccos{iz/i). The frequency of orbital plane precession 
in a triaxial potential, z/n, is ^ Vp^fi for low-i! orbits 
(pyramids and saucers) and even lower for other orbits 
(Figure [5|) . The corresponding timescale may be written 
as 

^ A^. P ,^^. 

in,triax ;C t.t r - (9') 



Comparison with the VRR timescale (|95p shows that 

iiifl.u/to.triax ^ ^fWt. For the Milky Way, these two 
timescales are roughly equal at a ~ 0.5 pc (Figure fT5|) . 
For sufficiently large IS! the regular precession due to tri- 
axial torques goes on faster than the relaxation, so the 
coherence time for VRR is now defined by orbit preces- 
sion, and the relaxation time itself becomes even longer. 
On the other hand, for small enough N the VRR destroys 
orientation of orbital planes before they are substantially 
affected by triaxial torques. It seems that VRR in triax- 
ial (or even axisymmetric) systems can be suppressed by 
regular orbit precession; we defer the detailed analysis of 
relaxation for a future study. 

So far we have considered the torques arising under 
RR as being independent of the torques due to the elon- 
gated star cluster. Suppose instead that we identify the 

"^ If the BH is spinning, pr ecession due to the Lensc-Thirring 
effect also destroys coherence IMerritt et al.|[201CII) . 



' N torques that drive RR with the torques due to the 
triaxial distortion. The justification is as follows: Dur- 
ing the coherent RR phase, the gravitational potential 
from N orbit-averaged stars can be represented in terms 
of a multipole expansion. If the lowest-order nonspheri- 
cal terms in that expansion happen to coincide with the 
potential generated by a uniform-density triaxial cluster, 
the behavior of orbits in the coherent RR regime would 
be identical to what was derived above for orbits in a tri- 
axial nucleus. We stress that this is a contrived model; in 
general, an expansion of the orbit-averaged potential of 
N stars will contain nonzero dipole, octupole etc. terms 
that depend in some complicated way on radius. Nev- 
ertheless the comparison seems worth making since (as 
we argue below) there is one important feature of the 
motion that should depend only weakly on the details of 
the potential decomposition. 

Equating the torques due to RR 



Tt 



RR 



N- 



Gr 



with those due to a triaxial cluster. 



Tt 



GNr 



triax 



our ansatz becomes 



7V-1/2. 



(98) 



(99) 



(100) 



As shown above (SJ7]), GR sets a lower limit to the 
angular momentum of a pyramid orbit (equation I74p : 



f . 



e a 



M{a) 



N- 



(101) 



the third term comes from setting (.q sa £„iax ~ >/e, the 
maximum value for a pyramid orbit, while the fourth 
term uses our ansatz (jlOOp and the definition (|65bl) of >c. 
Expressed in terms of eccentricity. 



1 



(^)' 



M, 



1 



N{a) 



(102) 



There is another way to motivate this result that does 
not depend on a detailed knowledge of the behavior of 
pyramid orbits. If we require that the GR precession 
time: 

_1 ^ _1 



(eq. 



'GR 



be shorter than the time 



dt 
for torques to change (. by of order itself, then 

»'Schw Af, 



< 



7V>f : 



M{a) 



N 



(103) 



(104) 



(105) 



as above. In other words, when I < ^min, GR precession 



is so rapid that the \/N torques are unable to change the 
angular momentum significantly over one precessional 
period. 

In order for this limiting angular momentum to be rele- 
vant to RR, the timescale for changes in the background 
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potential should be long compared with the time over 
which an orbit with i « ^min appreciably changes its an- 
gular momentum. As just shown, the latter timescale 
is 



tGR = i^ok ~ -"p' ~ ^Nvp\ (106) 



The former timescale is the coherence time for VRR, 
equation (|94l) : 

M. P 



The condition to ^ icR is then 
M. P 



N 



N ^ 



or 



a r- M, 
r-Schw m 



(107) 



(108) 



(109) 



Applying this to the center of the Milky Way, the condi- 
tion laecomes 

-^^/N(a)->li9 (110) 

mpc 

which is likely to be satisfied beyond a few mpc from 
SgrA*. 

On timescales longer than ^ icoh, the torques driv- 
ing RR will change direction. This is roughly equiva- 
lent in our simple model to changing the orientation of 
the triaxial ellipsoid, or to changing £o £^t fixed L Such 
changes might induce an orbit to evolve to values of £ 
lower than i'minj by advancing down the narrow "neck" 
in the lower left portion of Figure [T31 However such evo- 
lution would be disfavored, for two reasons: (1) it would 
require a series of correlated changes in the background 
potential, increasingly so as I became small; (2) as i de- 
creased and j^GR increased, changes in the background 
potential would occur on timescales progressively longer 
than the GR precession time, and adiabatic invariance 
would tend to preserve t (§4.4). These predictions can 
in principle be tested via direct A^-body integration of 
sma h-A^ systems i ncluding post-Newtonian accelerations 
(e.g lMerritt et al.llToiO) . 

A lower limit to the angular momentum for orbits near 
a massive BH could have important implications for the 
rate of gravitational wave events due to extrem e-mass- 
ratio inspirals, or EMRIs (jHils fc Bendeii I1995D . The 
critical eccentricity at which the orbital evolution of a 
IOM0 compact object begins to be dominated by gravi- 
tational wave emission is 



1 — gemri ~ 10 



U 



lO^yr 



-2/3 



\ IO^Mq 



4/3 



(111) 



with tr the relaxation time (e.g. lAmaro-Seoane et al.l 
120071 eq. (6)). By comparison, equation (|102p . after sub- 
stitution of N{< a) — A''o(a/mpc) implies 



1-en 



2x10 



-if M, 



V IO^Mq 



100 



-1 / 



\ lOmpc 



(112) 



9. ESTIMATES FOR REAL GALAXIES 

In this section we estimate the fraction and lifetime 
of pyramid orbits to be expected in the nuclei of real 
galaxies. 

We restrict calculations to the case of "maximal triax- 
iality," £& = ec/2, although we leave the amplitudes of 
eb, €c free parameters. We also limit the discussion to or- 
bits within the BH influence radius, r < ?'infl, where our 
analysis is valid and where orbits are typically regulai[j. 
Beyond ^ rinfi, centrophilic (mostly chaotic) orbits still 
exist and c ould dominate, e.g., th e rate of feeding of a 
central BH (jMerritt fc Poonll2004D . 

The first set of parameters is chosen to describe the 
center of the Milky Way. The BH mass is set to 
M, = 4 X IQS Mq iJGheTet all 120081 : iCillessen et al.l 
l2009al |a) . The density of the spherically symmetric stel- 
lar cusp is taken to b e ps = 1.5 x lO ^Mfr^pc"^ (r/1 pc)~ '^ 
(|Schodelet al.ll2ml . with 7 = 1.5 (jSchodel et al.ll2008fh 
the corresponding BH influence radius is n^d ~ 3 pc0 
The triaxial component of the potential is highly un- 
certain; one source wo uld be the nuclear bar with den- 
sity p t ~ 15OM0pc~^ (jRodriguez-Fernandez fc CombesI 
l2008f ). yielding a triaxiality coefficient at r = 1 pc of 
Cc ~ 10"^. We also considered a larger value, €c — 10^^, 
which may be justified by some kind of asymmetry on 
spatial scales closer to rjnfl than the bar. In this model, 
the precession time due to the spherical component of the 
potential, 27r/(3fp) for a circular orbit, is independent of 
radius and equals '-^ 1.7 x 10^ yr; the two-body relax- 
ation time is also constant (5 x 10^ yr), and timescales 
for scalar and vector resonant relaxation are given by 
equation ([M]) . &5\i with stellar mass m = 1 Mq. 

The second set of parameters is intended to describe 
the case of galaxies with more massive BHs, using the 
so-called M, — a relation in the form 



A/, w 1.7 X 10* Mq 



200 km s"^ 



4.86 



(113) 



(jFerrarese fc FordI 120051 ). Combined with the definition 
of Tinfl « GM./cr^, we get 



Tinfl W 13 pc 



M, 



108 Mq 



0.59 



(114) 



The two-body relaxation time evaluated at rjnfl (as- 
suming a mean-square stellar mass rrii, = 1 Mq and a 
Coulomb logarithm In A = 15) is 



<26r(n„fl)~2.1x 10^^ yr 



200 km s" 



7.5 



^9.6 X 10' 



yr 



M, 



IO^Mq 



1.54 



(115a) 



(115b) 



* Excepting for the effects of GR, which as noted above may in- 
troduce chaotic behavior even for orbit-averaged parameters. The 
chaos that sets in at r > rjufj (§[6ll arises from the couphng of the 
orbit-averaged and radial motions. 

^ These cusp parameters correspond to an inward extrapola- 
tion of the de nsity o bserved at r > 1 pc. Recent obs ervations 
HBuchholz et aI.ll20C)9l : IDo et al.|[200g : IBartko et aI|[20Tg ) reveal a 
"hole" in the density of evolved stars inside ~ 0.5 pc, implying a 
possibly much lower density for the spherical component near the 
BH. 



Orbital structure of triaxial nuclei 



19 




Fig. 15. — Various timescales for Milky Way center. Green and 
blue solid curves are the pyramid draining times (more exact cal- 



culation than equation I117|l for Ec 



10" 



and 10 ; dotted green 



and blue curves denote typical orbital plane precession timescales 
u^ for low-£ orbits (£^ ^ ec{a)); Newtonian and relativistic pre- 
cession times are given for circular orbits; other timescales are 
marked on the plot. Vector resonant relaxation is suppressed when 
i'RR,v f^ i^Q (marked by conversion of line into dashed). 

(jMerritt et al.ll2007t ). 

We first estimate the radius rcrit that separates the 
empty {q < 1) and full loss cone regimes. As noted in 
the previous section, GR precession prevents a pyramid 
orbit from reaching arbitrarily low angular momenta; the 
radius beyond which capture becomes possible is roughly 
Tcrit. Using equation (f56| with W = (15ec)^ (the maxi- 
mum value for pyramids) and equations (|14a|) . (|44|) . the 
condition q = 1 translates to 



1 



An psr^ /rcrit 



3-7 M. 



ro 



3-7 



2a' 



5 Tree (rcrit) 

3(2-7) V0^schw/r, 



crit 



If we take tq to be rjnfl, and ps and a as the density and 
velocity dispersion at this radius, we obtain 



2/7 



rcrit 

r-infl 



0.5^ 



c ec(ro) 



(116) 



The radius rcrit typically lies in the range (0.2 — 0.7)rinfl, 
weakly dependent on the parameters. Since regular pyra- 
mid orbits exist only for a < r-mfi (Figure ITOl) , there is 
evidently a fairly narrow range of radii for which cap- 
ture of stars from pyramid orbits is possiblcl3 However 
pyramid- lik e, ccntrophilic can exist at much larger radii 
(|Poon fc Mcrritt 200li\ 

Next we make a rough estimate of the pyra mid dr ain- 
ing time at a > rcrit, using the expression (j60bp for 
the flux in the full loss cone regime, J" ~ p,/{Pvp); p, 

^^ This is also roughly the radial range from which e xtreme- 
mass -ratio inspiral events are believed to originate; (e.g. Ilvanovl 
[2001) . 



(the fraction of phase space occupied by the loss cone 
is given by equation ([51]) with £ ■ -^ " ^^ ~ ^ 
(equation [SU]) : 






idr 



biT 



Tv„ 



' 10 yr X 



39 (GM.)3/2 
ec(a) f M, 



,5/2 



e V 108 Mr, 



eda) 



-3/2 



1 pc 



(117) 

5/2 



A more exact calculation of idrain(a) for the Milky Way, 
based on numerical analysis of properties of orbits sam- 
pled from the entire phase space, is shown in Figure [T5l 
Finally, we estimate the total capture rate for all pyra- 



mids inside ri„fl, using t 



£W 



timescale and applying (|1131 1114[ I117P : 
ec(ri„fi) f M, 



^drain(rinfl) as a typical 



t 



pyr 



6 X 10 yr x 



The capture rate from pyramids is then 

ecirinfi)M, 



-0.025 



(118) 



M, 



pyr 



t 



1.6xlO"'*M0 yr^i 



e 



pyr 



M, 



1.025 



108 Mq 

(119) 

For the Milky Way we find ~ 4 x 10 ^ Mqjt ^ for Ec = 
10-3 and ~ lO-^'Mgyr-i for e^ = 10"^. 

This capture rate should be compared with tha t due to 
two-b ody relaxation, which is estimated to be (jMerrittl 
[200l 



M, 



2br 



0.1 



Ml. 

t2br 



lO-'^Afoyr^i x 



M. 



108 M, 







-0.54 



(120) 

Thus even for a Milky Way-sized galaxy, the capture 
rate of pyramids could be comparable with or greater 
than that due to two-body relaxation. For more massive 
galaxies this inequality becomes even stronger. However, 
this is only the initial capture rate - after ~ ipyr, all stars 
on pyramid orbits would have been consumed, at least 
in the absence of other mechanisms for repopulating the 
small-.^ parts of phase space (not necessarily £ ^ i,, but 
the much broader region £ < y^e^ from which draining is 
effective) . 

In the most luminous galaxies, like M87, standard 
mechanisms for relaxation are expected to be ineffec- 
tive even over Gyr timescales and pyramid orbits once 
depleted are likely to stay depleted. Setting Ec — 0.1, 
e = 3 and A/. = 4 X IO'^Mq gives for M87 ipyr « 5 
Gyr and Mpyr « 4 x IO^Mq. This could be an effective 
mechanism for creating a low-density core at the centers 
of giant elliptical galaxies. 

10. CONCLUSIONS 

We discussed the character of orbits within the radius 
of influence rinfl of a supermassive BH at the center of 
a triaxial star cluster. The motion can be described as 
a perturbation of Keplerian motion; we derive the orbit- 
averaged equations and explore their solutions both ana- 
lytically (when the triaxiality is small) and numerically. 
Orbits are found to be mainly regular in this region. 
There exist three families of tube orbits; a fourth or- 
bital family, the pyramids, can be described as eccentric 
Keplerian ellipses that librate in two directions about 
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the short axis of the triaxial figure. At the "corners" of 
the pyramid, the angular momentum reaches zero, which 
means that stars on these orbits can be captured by the 
BH. We derive expressions for the rate at which stars on 
pyramid orbits would be lost to the BH; there are many 
similarities with the more standard case of diffusional 
loss cone refilling, but also some important differences, 
due to the fact that the approach to the loss cone is 
deterministic for the pyramids, rather than statistical. 
The inclusion of general relativistic precession is shown 
to impose a lower bound on the angular momentum. We 
argue that a similar lower bound should apply to orbital 
evolution in the case that the torques are due to resonant 
relaxation. The rate of consumption of stars from pyra- 
mid orbits is likely to be substantially greater than the 



rate due to two-body relaxation in the most luminous 
galaxies, although in the absence of mechanisms for or- 
bital repopulation, these high consumption rates would 
only be maintained until such a time as the pyramid or- 
bits have been drained; however the latter time can be 
measured in billions of years. 
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